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Abstract 

We solve numerically the equations of nonlinear fluctuating hydrodynamics (NFH) for the su- 
percooled liquid. The time correlation of the density fluctuations in equilibrium obtained here 
shows quantitative agreement with molecular dynamics (MD) simulation data. We demonstrate 
numerically that the 1/p nonlinearity in the NFH equations of motion is essential in restoring 
the ergodic behavior in the liquid. Under nonequilibrium conditions the time correlation func- 
tions relax in a manner similar to that observed in the molecular dynamics simulations in binary 
mixtures. The waiting time t w dependence of the non-equilibrium response function follows a Modi- 
fied Kohlrausch- Williams- Watts(MKWW) form similar to the behavior seen in dielectric relaxation 
data. 
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The conserved densities of mass, momentum, and energy constitute the simplest set of 
slow modes characteristic of the isotropic liquid. The microscopic balance equations for the 
respective conservation laws contain terms with widely different characteristic time scales 
of variation respectively corresponding to the various degrees of freedom of the complex 
system. The nonlinear fluctuating hydrodynamics (NFH) describes the dynamics of these 
slow modes with nonlinear differential equations having regular and stochastic parts. The 
regular parts involve nonlinear coupling of slow modes while the random parts represent 
noise which can be linear [l|, [3, [J or multiplicative^, 4|. The most widely studied theoretical 
model for the slow dynamics in a supercooled liquid approaching vitrification follows from 
;hese equations of NFH and is termed as the self-consistent mode coupling theory (MCT) 
5), [f| . In a strongly interacting dense liquid the coupling of density fluctuations produces 
the dominant effect on dynamics. The MCT involves a nonlinear feedback mechanism [3] 
of density fluctuations producing strong enhancement of the viscosity of the supercooled 
liquid. In its simplest version the MCT predicts that above a critical density the long 
time limit of the time correlation C(t) of density fluctuations is nonzero. This signifies 
an ergodic-nonergodic transition(ENE) in the liquid and is a precursor to the liquid-glass 
transition. The predicted dynamics involves several different regimes of relaxation and has 
been widely used in fitting experimental data on different liquids. However, the simple MCT 
approach is known to exaggerate the tendency of the dynamics towards slowing down, and 
becomes quantitatively inaccurate in the vicinity of the predicted transition, which is never 
observed in practice. The perturbation expansion for the renormalized transport coefficient 
in the MCT, though systematic, is in terms of a dimensionless parameter which is not small. 
Furthermore, it has also been shownQ, Q| non perturbatively that the 1/p nonlinearities in 
the NFH equations remove the sharp ENE transition predicted in the simplified theory. We 
report here the study of the slow dynamics of a dense monatomic Lennard- Jones liquid by 
numerically solving the stochastic equations of NFH. Our nonperturbative calculation shows 
good agreement with the computer simulation results of the same system in equilibrium. We 
also report here results on the dynamics of the density fluctuations under nonequilibrium 
conditions. 

We consider the equations of NFH for an isotropic liquid in its simplest form for the mass 
density p and momentum density g[l| 
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The correlations of the gaussian noise 9i are related to the bare damping matrix L^- 
(e i (x,t)e j (x't')) = 2k B TL ij 5(t-t')5(x-x / ). For an isotropic liquid, L tj = (Co + W3)<%V 2 + 
r/oViVj where Co and r/o respectively denote is the bare bulk and shear viscosities. The 
stationary solution to the Fokker-Planck equation corresponding to the generalized Langevin 
eqn. is obtained as exp{— /3F[p, g}} with /3 = 1/ksT is the Boltzmann factor. The coarse 
grained free energy functional is obtained as F[p, g] = Fjc[p,g] + Fjj. The kinetic part is 
dependent on the momentum density F K = J dxg 2 / (2p) and the so called potential part is 
given by F v = F id + F int . The ideal gas contribution is F id = f drp(r)[ln(p(r)/p ) ~ !]■ The 



interaction part F int up to quadratic order in density fluctuations [lOj is obtained as 



pFtot = ~2~2 J drdr ' c ( r - r')Sp(v)5p(r') , (3) 

n 

where c(r) is the two point Ornstein-Zernike direct correlation functional and m is the mass 
of the particles. For the glassy dynamics we focus on the coupling of slowly decaying density 
fluctuations present in the pressure functional, represented by the third term on the LHS of 
eqn. (j2J). With the above choice of Fjj, the nonlinear contribution in this term reduces to 
pVif(r,t) with the convolution f(r,t) = rrT 1 J drc(r — r')5p(r ,t). 

We consider here a classical system of N particles, each of mass m interacting with a 
Lennard- Jones potential u(r) = 4e[(cr/r) 12 — (cr/r) 6 ]. In addition to the scale a of the 
interacting potential there is another length h of the lattice grid on which p and g are 
computed. We choose a/h to be non integer ( = 4.6 in the present calculation) to avoid 
crystallization. Time is scaled with the LJ unit of r = {ma 2 /e)a and length with h. The 
thermodynamic state of the fluid is described in terms of the reduced density n* = no<j 3 
and the reduced T* = (fcsT)/e. For numerical solution the conserved densities are scaled 
to dimensionless forms: n(r) = [/i 3 m _1 ]p(r), and j(r) = [/i 3 (me)^%(r). The speed of 
sound Co is given by, Cq = ksT '/ \mS(Q)). We start with an initial distribution of the 
fluctuating variables n(r) and j(r) over a set of points 20 3 on a cubic lattice. The nonlocal 
integral f(r,t) is evaluated as a sum of contributions from the successive shells, f(r,t) = 
h^Yli c (Ri)Ylia^ n (R'i it), where Rf for a = 1, ...rrii respectively denote radii vectors of 



the rrii lattice points in the zth spherical shell of radius R{. The 1/p nonlinearity in the 
dissipative term of the momentum equation is computed by replacing the density field in the 
denominator with the p(x) averaged over a length scale close to o around the corresponding 
point r. We ignore the convective nonlinearity in the present calculation and focus on the 
role of the pressure nonlinearity in producing the slow dynamics. 

A major hurdle encountered in the numerical scheme used here arises from an instability 
which occurs as n(x, t) gets negative at certain grid points. To avoid this situation, we 
redefine n(x, t) on the grid at each step of the numerical integration with a coarse graining 
scheme. In devising the latter we make use of the following physical interpretation of the 
definition of p(x, t) of the density field : the integral f AV <ixp(x, t) represents the total mass 
in an elementary volume AV of the system. At each time step of the numerical integration, 
the positivity of the field n(x) over the whole grid is checked. If it turns negative at a point, 
we reduce n(x) at some or all of the neighboring sites by taking equal contributions from 
each and add the sum total to the original site. It is also ensured that the density at none 
of the neighboring sites becomes negative as a result of this redistribution. The sum of 
the densities at the original and the contributing sites remains unaltered and hence global 
conservation is maintained. If the above redistribution involving contributions only from 
the nearest neighbor sites is insufficient to make n(x) positive everywhere, we include the 
next nearest neighbors in the redistribution and so on. In reality however we hardly need to 
include beyond the second shell of neighbors surrounding the original site. With the density 
instability being corrected with this coarse graining procedure, the numerical algorithm can 
be run up to much larger times than in earlier works HI] . The arbitrary regularization of 
the strength of the noise llj can also be avoided, and the fluctuation dissipation relation 
respected. 

The equal-time correlation of density fluctuations for the N particle system is given by 
S(k,t) = N^ 1 < 5n(k,t)Sn(—k,t) > where the angular brackets refer to an average over 
the noise. We first consider the system as it evolves at T* = 2.0 and = .97 under the 
influence of thermal noise, starting from an initial state in which all fluctuations are set 
to zero. Time translational invariance is reached as the system equilibrates and S(k, t) 
approaches the corresponding static structure factor S(k). The equilibrated S(k) vs. k plot 
( for large t w ) is displayed in fig. 1. S(k,t) obtained with equations of motion linear in 
fluctuations is also displayed. The peak position (q m ) and amplitude of the S(k) obtained ( 
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FIG. 1: S(k) vs ka at T* = 2 and Hq = 0.97 for linear (solid) and nonlinear (dashed) dynamics; as 
computed from input c(r) (dotted). Inset: S(k,t) vs t for ka = 6.75 (solid) and k = 7.05 (dashed) 
displaying equilibration with time. 



using the Ornstein-Zernike relation) from the input direct correlation function c(r) are well 
reproduced. Other features of S(k) are partly lost due to the relatively crude grid size used 
in our numerical solution. 

Next we focus on the dynamic correlation function C(t + t w ,t w ) defined in the normalized 
form C(t+t w , t w ) = < 5n(t + t w )5n(t w ) >/< 5n(t w )5n(t w ) >. For large t w time translational 
invariance holds, making C(t+t w , t w ) = C(t). The decay of C(q m , t) is compared in fig. 2 with 
the corresponding molecular dynamics simulation results 12|] for the equilibrated systems at 
T = .6 for two densities Uq = 1.10 and = 1.06. The input c(r ) corresponds to the purely 



repulsive part of the Lennard- Jones potential following Ref. [12[]. The bare transport coeffi- 
cients which determine the noise correlations are chosen such that the short time dynamics 
agrees with computer simulation results. C(t) obtained by solving the stochastic equations 
linearized in the fluctuations decays very fast. In comparison considerable slowing down of 
the decay of C(t) occurs on solving the full NFH equations. The static correlation function 
S(k) however ( see inset of fig. 3) shows hardly any difference between the two cases. At 
higher densities the mean-free path of the fluid particles gets smaller and approaches the 
atomic length scale. As a result, the validity of Generalized hydrodynamic equations at 
short length scales ( corresponding to wave vector g~g m ) improves with increasing density. 
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This trend is clearly seen our results displayed in fig. 2 for C(t). 

The ENE transition of simple MCT is driven by the nonlinear couplings of density fluc- 
tuations in the pressure term (3rd term on LHS) of the generalized Navier-Stokes equation 
On the other hand the 1/p nonlinearity crucial for the absence of the ENE transition 
is in the dissipative term of the same equation, i.e., 4th term on LHS of ([2]). We therefore 
consider two cases here to test the role of the relevant nonlinearities from a non perturba- 
tive approach. In Case A the 1/p nonlinearity in ([2]) is replaced with 1/po while keeping 
the density nonlinearity in the pressure term. In Case B, the complete model with both 
nonlinearities is considered. The results for C(q m ,t) at T = .6 and Uq = 1.10 corresponding 
to the Cases A and B respectively are shown in fig 3. We extend the numerical solution to 
the longest possible time scale ( > 10 3 in Lennard- Jones units) which is about four orders 
of magnitude beyond the microscopic time scales. The decay of the dynamic correlation is 
markedly different in the two cases and agrees with the previous theoretical results on the 
role of 1 j p nonlinearity. 

To study the structural relaxation in the nonequilibrium state we consider the time 
evolution of the Lennard- Jones liquid following an instantaneous quench from T* = 2.0 
and Uq = 0.97 along the isobaric line to Tj = 0.4 and Uq = 1.12. We compute from the 
solution of the NFH equations the two-time density correlation function C(t w + t,t w ) at 
q = q m for different waiting times t w = 50, 100, 200, 500, 700 and 1000. For small values 
of t time translational invariance holds and C(t w + t, i w ) depends only on t. On the other 
hand at large t, the correlation function depends on both t and t w . Following the mean-field- 
theoretical results 13] and also experimental datajjjl fits on spin glasses, we fit this long time 
part of the density correlation function with the form C ag [h(t + t w )/h(t w )}, where h(t) is a 
monotonously ascending function of its argument. In the C (i w + 1, £ w ) vs. [h{t + t w ) / h(t w )] a 
plot, the parameter a is tuned to obtain a collapse of all the curves. The results displayed 
in the fig. 4 obtain h(t) ~ [log(t)] a with the best fit value of the parameter a = .81. This is 
comparable to the corresponding value a = .88 obtained in molecular dynamics simulations 
{15I ] of binary Lennard- Jones mixtures. 



In a recent work, Lunkenheimer et. al. 



16| studied over a range of frequency u the 



dielectric response function Xw(^w) at temperature T < T g . The system falls out of equilib- 
rium over laboratory time scales at the calorimetric glass transition temperature T g . The 
aging time (= t w in the present notation) dependence of x follows a modified Kohlrausch- 
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Williams- Watts (MKWW) function f(t w ) = exp[(t w /r(t w )) /3 ]. The relaxation time r(t w ) and 
the stretching exponent f3 are identical for all frequencies. The limiting value r(t w — > oo) is 
close to the a-relaxation time r a extrapolated to corresponding temperature T < T 9 [ly]. In 
the present work T ~ T c (> T g ) and in this case the system in fact equilibrates. We study 
here the function Xuj(t w ) = LuC(uj,t w ) which in equilibrium would reduce to the correspond- 
ing response function. C(u,t w ) is obtained approximately (i.e. ignoring FDT violations) 
from the frequency transform of C(t + t w ,t w ) with respect to t. The data is fitted to the 
form : 

xM = [x s J-x e J}f(t w ) + x e J, (4) 

where x^ an d xt? respectively refer to the initial and final values of Xw F° r the relaxation 



time in / we usell7(] r(t w ) = (r st — r eq ) f(t w ) + r eq with the normalized function f(t) 
2^/[l + exp{2t/r(t)}] /3 . Fig. 5 shows that x^'s at different frequencies scale onto a single 
master curve for (3 = 0.68. The r(£ w ) decreases sharply with t w initially and becomes almost 
constant at r cq . At large t w the relaxation follows a stretched exponential form having the 
a-relaxation time r eq and exponent j3 at the corresponding temperature. This shown in the 
inset of fig. 5. 




FIG. 2: C(q,t) vs t/ro at T* = .6 for densities = 1.10 (solid) and 1.06 (dotted). Corresponding 
MD simulation data[12] shown respectively with (dashed) and (dot-dashed) curves. Inset : short 
time part of C(q,t) with (solid) and without (dashed) nonlinear coupling of modes; corresponding 
MD simulation result [12] (dark circles). 
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FIG. 3: C(q,t) vs t/rg at T* = .6 and Uq = 1.10 for case A (filled circles) and case B (filled 
triangles). The solid lines are the best fit curves to the corresponding data. Inset : S(k) vs ka for 
the two cases respectively with solid and dashed curves. 
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FIG. 4: Normalized correlation function C(t w + i,t w ) vs. [h(t + t w )/h(t w )] a for different waiting 
times i w / T o shown in the inset . The data at large t scales onto a single master curve with 
h(t) ~ log(t) and a = .81. 



We have shown that the direct numerical solutions of the NFH equations provide a 
reliable way of studying the dynamics of fluctuations in a dense liquid in the vicinity of the 
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FIG. 5: Scaling of Xw(*w) f° r frequencies ujtq = 10 _2 (circle), 10 _1 (diamond), 10° (triangle), and 
10 1 (star), solid line is a fit to MKWW form. Lower inset : r(t w ) vs iwAo- Upper inset : C(q m ,t) 
vs i/ro at the final state T* = 0.4 and = 1.10 (circles); solid line is fit to a KWW form with 
exponent (3 = 0.68 and relaxation time r eg (see text). 



avoided ergodic nonergodic transition. While the method is not appreciably more efficient 
than molecular dynamics from a computational standpoint, it provides an interesting way of 
investigating theoretical assumptions that can be made in the analytical treatment of these 
equations. In particular, our work clearly shows the role of the 1/p nonlinearity present 
in the dissipative term of eqn. ([2]) in restoring ergodicity [lj. The present method can be 
easily extended to a larger set of hydrodynamic variables, which would permit a description 
of the dynamics in binary mixtures. For the nonequilibrium states, the dynamics with the 
one loop mode coupling theory has been formulated for the spherical p-spin glass model [3]. 
For the supercooled liquid a similar analysis even at the one loop level is still lacking. Our 
numerical approach is by essence non perturbative, and it will be interesting to compare it to 
analytical perturbative approaches to the same equations under nonequilibrium conditions. 
The CEFIPRA is acknowledged for financial support under Indo-French research project 



S 



2604-2. 



[1] S. P. Das and G. F. Mazenko, Phys. Rev. A, 34, 2265 (1986). 

[2] A. Andreanov, G. Biroli, and A. Lefevre, J. Stat. Mech, PO7008 (2006). 

[3] T. N. Nishino and H. Hayakawa, Cond-mat. 0803.1797v2 

[4] B. Kim and K. Kawasaki, J. Stat. Mech, P02004 (2008). 

[5] S. P. Das, Rev. of Mod. Phys. Rev. 76, 785 (2004). 

[6] D. R. Reichman and P. Charbonneau, J. Stat. Mech., P05013 (2005). 

[7] U. Bengtzelius, W. Gotze and A. Sjolander, 1984, J. Phys. C 17, 5915. 

[8] S. P. Das and G. F. Mazenko, Cond-mat 0801.1727. 

[9] J-P. Hansen and I.R. . McDonald, Theory of Simple Liquids, 3rd ed. Academic, London (2006). 

[10] T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B. 19, 2775 (1979). 

[11] L. M. Lust, O. T. Vails and C. Dasgupta, Phys. Rev E. 48, 1787 (1993); O. T. Vails and G. 

F. Mazenko, Phys. Rev A. 46, 7756 (1992). 

[12] J. J. Ullo and S. Yip, Phys. Rev. Lett 54, 1509 (1985) 

[13] J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and M. Mezard, Physica (Amsterdam) 226A, 

243 (1996). 

[14] H. Rieger, J. Phys. A 26, L615 (1993). 

[15] W. Kob and J. L. Barrat, Phys. Rev. Lett. 78, 4581 (1997). 

[16] P. Lunkenheimer, R. Wehn, U. Schneider, and A. Loidl, Phys. Rev. Lett. 95, 055702 (2005). 

[17] B. Sengupta and S. P. Das, Phys. Rev. E, 76, 061502 (2007). 

[18] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 173 (1993). 



9 



